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Abstract This paper explores the dynamic properties of the planar system of an ellipsoidal satellite in 
an equatorial orbit about an oblate primary. In particular, we investigate the conditions for which the 
satellite is bound in librational motion or when the satellite will circulate with respect to the primary. We 
find the existence of stable equilibrium points about which the satellite can librate, and explore both the 
linearized and non-linear dynamics around these points. Absolute bounds are placed on the phase space of 
the libration-orbit coupling through the use of zero-velocity curves that exist in the system. These zero- 
velocity curves are used to derive a sufficient condition for when the satellite's libration is bound to less 
than 90°. When this condition is not satisfied so that circulation of the satellite is possible, the initial 
conditions at zero libration angle are determined which lead to circulation of the satellite. Exact analytical 
conditions for circulation and the maximum libration angle are derived for the case of a small satellite in 
orbits of any eccentricity. 

Keywords Libration • Gravity Gradient ■ Binary Asteroids ■ Full Two-body Problem ■ Libration-Orbit 
Coupling 



1 Introduction 

The investigation of the translational-rotational coupling for a finite orbiting body, referred to in the 
literature as the Full Two Body Problem, has received renewed attention in recent years |2 H ll H[T0l Hlll6]. 
Much of this work has been motivated by interest in the dynamic evolution of binary asteroid systems, 
which comprise 16% of Near- Earth asteroids [T2] . Scheeres jT5] studied the stability of bodies resting on 
one-another which can lead to the formation of a binary asteroid system through rotational fission. Bellerose 
[51 looked at the dynamics and stability in the planar Full Two Body Problem from an energetic standpoint. 
Fahenstock 6 studied the problem by modeling the gravitational interaction with polyhedral models, which 
were applied to simulate the dynamics of the near- Earth binary asteroid 1999 KW4 [7]. 

Over time, many simplifications have been made to the Full Two Body Problem in order to make 
analytical progress in the study of specific applications. The most common simplification is to treat the 
finite bodies to second-order in their mass properties, which allows for the use of inertia matrices for 
modeling the attitude dynamics [9ll2l ll9l [3]. Analytical expressions have been constructed to fourth order 
|20| . however most analytical studies stop at second order. A second simplification comes by reducing the 
problem to the planar problem. If the second order simplification is made on mass properties, then the 
in orbit plane libration angle (commonly referred to as the pitch angle) dynamics are decoupled from 
the out-of- plane (roll/yaw) dynamics if they are initially quiescent [9j. This reduces the dimensionality 
of the problem from six to three degrees-of-freedom. The third major simplification that is often seen, 
particularly in the spacecraft community, is decoupling the translational and rotational dynamics [9l[2]. In 
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this case, it is assumed that the mass of the orbiting body is insignificant compared to the primary, so 
that the perturbation to the orbital dynamics is neghgible. The orbital dynamics are known from the given 
Keplerian motion, and these are used as inputs to the attitude dynamics. A notable work which explores 
the effects of the coupling for the case of spacecraft sized objects was carried out by Mohan IS . 

In this paper, we analyze the problem of a triaxial body in orbit about a spherical primary, or in an 
equatorial orbit about an oblate primary. These bodies are represented to second order with their moments 
of inertia, but all coupling between the translation and attitude motions for this case are preserved. This 
model can be used as a first order representation of common dynamic problems including binary asteroids 
(secondaries tend to be in near equatorial orbits of oblate primaries [H]), spacecraft orbiters, and planet- 
moon systems. This paper is closely related to those by Scheeres [19] and Bellerose [3], however we extend 
the foundation laid by those works. 

The main contribution of this paper is the study of the limits of librational motion for a given system. 
There is a small amount of literature dedicated to bounding librational motion. Auelmann 1 laid out the 
bounding conditions for an axially symmetric spacecraft in circular orbit, using the uncoupled attitude 
dynamics with a zero spin rate about the axis of symmetry. Pringle 15J extended this idea, analyzing the 
equilibrium cases and librational bounds for all spin states. 

This paper only looks at the planar libration, however we account for coupling, triaxial shapes, and are 
not limited to circular orbits, or indeed even Keplerian orbits due to the perturbations from the attitude 
coupling. The methodology used to bound the libration is initially similar to [T1 I15] . in that we use the 
energy to limit the amount of libration that is dynamically possible by computing zero-velocity curves 
for the system. We extend this, however, to cases when the energy is high enough so that circulation is 
possible. In this case, we show that although the energy zero-velocity curve does not bound the libration 
angle, certain initial conditions can lead to trajectories which have bounded maximum libration angles. In 
particular, for the case where the ellipsoid is very small compared to the other body (such as a spacecraft 
in Earth orbit), we can derive analytical expressions which predict the maximum libration amplitude for 
an orbit of any eccentricity. For cases with significant coupling, we present a method of analysis which 
determines the structure of bound and unbound trajectories. 

The paper is organized as follows. Section [2] reviews the mathematical description of the system in 
question as originally derived by Scheeres [19 . Section [s] discusses the location and properties of equilibria 
in the system. Section |4] presents a sufficient condition for bounded librational motion. The fact that this 
is only a sufficient condition is shown in Section [5] Finally, Section |6] presents the main results of the paper 
which describe conditions for unbound libration. 



2 Physical System Description 

In general, we begin with the planar full two-body problem where both bodies are finite, but constrained 
to have their individual rotation poles perpendicular to the mutual orbit plane. This situation is illustrated 
in Fig. [T] with the angles of interest and radius defined. 
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Fig. 1 Definition of tlie system angles and relative radial vector. 
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In this section we present the key mathematical relationships used to describe this system, which are 
simplifications of the general relationships derived in [19]. Further details of the derivations are available 
in the Appendix. 



2.1 Equations of Motion 

Scheeres [1_QJ showed that in this case, the potential energy using a second order expansion in the moments 
of inertia is, 



V{r,(j>i,cj>2) 



GMiAh 



1 + 
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2^ 
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where I,^ is the mass normalized (signified by the prime) inertia dyad of body i, and the subscripts on the 
non-bold versions indicate a principle moment of inertia. G is the gravitational constant, Mi the mass of 
body i, Tr() indicates the trace of the dyad. 

In this paper, we will normalize the lengths with respect to the maximum ellipsoid semi-axis, a; the 
time will be normalized by the mean motion of the system at this distance, n = \J G{Mi + M2)/cfi (units 
of 1/s); and the mass will be normalized by the ellipsoid mass, Af2. Note that using these normalization 
factors, the normalized value of /i = G{Mi + M2) is 1. All equations from this point on will be written in 
normalized units. 

First, let us define the reduced mass as. 



the mass fraction is then. 



M1M2 
Ml + M2 

Ml 



Ml +M2 

which is just the reduced mass in normalized units. Next, we define the general moment of inertia as, 

Iz = /2^ + vr^ 



(2) 
(3) 
(4) 



Note that the general moment of inertia is a function of r. The normalized moment of intertias are computed 
as 

M2a2 a2 



(5) 



Mia2 a2 

For a system with an oblate primary, we have that Ii^ 
becomes in normalized units, 



(6) 



Iiy = Is, where Is < Ii^, so the potential 
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In the case of a spherical primary, which means that 
simplifies to. 



7i^ , and therefore the potential 
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Note that in both the spherical and oblate primary case, the spherical symmetry causes the potential 
energy to no longer be dependent on the primary's orientation, as represented in the planar problem by 0i. 
Starting from Scheeres' |19j expression, in this case the kinetic energy is. 
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where 61 = 4>i + 6. Note that the dot indicates derivatives with respect to unit-less time. 
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2.2 Integrals of Motion for the System 



The system described in Section [27T] has three integrals of motion that can be used to simplify the problem: 
the Jacobi constant (or total energy), the total angular momentum {Ktot), and the inertial angular velocity 
of the primary (^i). That these quantities are integrals is shown in the Appendix. 
Using the integrals, the kinetic energy from Eq. ^ can be rewritten as 



(10) 



where Ti = {Mi/(2M2))Ii.0\ is the kinetic energy of the primary, which is constant. Likewise, we can 
define the free angular momentum from Eq. (Ill to be. 



K = Kt. 



Ki = he + I2 



(11) 



where Ki = {Mi/M2)Ii^6i is the angular momentum of the primary, which is constant. This relationship 
between the free angular momentum and the orbit angular velocity will allow us to eliminate 9 from the 
system. Solving Eq. (Ill for the angular velocity gives. 



(12) 



Substituting Eq. ( 12 ) into Eq. ( 10 1 allows us to write the kinetic energy as, 

nn nn ^^'^ 1-2 1 ^2 I^v"^ <i>^ 

T = Ti H h -i^r H — — 

2 ^ 2 ^2 h 

so that the total energy of the system can be written as, 

r2 



2 12 
2^vr 02 



P ^ , 1-^ , 1 -2 , ll2^i^r 
or alternatively the free energy of the system can be written as, 

^2 



+ V{r,<l,2) 



212 



(13) 



(14) 



(15) 



It is interesting to note that the description of the free energy of the system does not require any 
knowledge about the primary spin state as neither 01 or 01 appear anywhere in Eq. (151. 



2.3 Dynamic System Analysis 

In order to look at the stability of any relative states for the planar system, we derive the dynamic matrix 
corresponding to the equations of motion 

A=^^ (16) 

where q = [r 6 cj>i 02] ■ 

First, however, we note that we can reduce the order of the system by recalling from Section [2. 2| that 
01 is ignorable, and that through the conservation of angular momentum we can remove 6 through the use 



of Eq. (12 1. Therefore the equations of motion for the 2 degree-of-freedom system are, 

^ ^ {K-l2A2fr IdV^ 
dr 



V Or 



(17) 



1 dV 2r{K-l2M) 



1 + ^^^+ ^ r (18) 

i2^ ) "02 rlz 

This system is depicted from an ellipsoid fixed frame in Fig. ([2]). 

This frame rotates with the orbit, so that if there is no libration of the secondary, the primary location 
will be fixed. When the secondary is librating, the primary will appear to move relative to the secondary 
due to the libration, as well as changing the separation distance as energy is traded between the orbit and 
the secondary libration. 
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Fig. 2 Definition of the secondary fixed relative frame. 



The dynamic matrix is now a 4x4 matrix with the form, 






dr 



1 
1 

#■ 



^ d(f)2 

d4>2 d<p2 d4>2 
dr d4>2 dr 



(19) 



where the state vector has been organized as, x = [r (t>2 r 4>2\. This matrix will be used in the following 
sections to analyze the stability of the system around any relative equilibrium points that may exist. 

3 Equilibrium Conditions 



In this section, we determine the relative equilibria between the two bodies. This means that these equilibria 



are places where the relative dynamics given by Eqs. (17) and (18) are stationary. However, the primary 



is free to rotate at a constant rate, and the secondary (ellipsoidal) body will orbit at a constant radius 
and orbit rate {8). As in Section [2] many of the relationships given here are simplifications of the general 
relationships derived in 19]. Further details are available in the Appendix. 



3.1 Relative Equilibrium Point Locations 



The relative equilibrium locations are found by solving the following polynomial for r. 



6 ^ 5^ 

r 5-r + 



2L 



+ 



4 + ^ (C± +7,^.-7. 



(20) 







(21) 



where 

' -2l2^+l2y +72, if 02 =0, TT 
- 272„+72, if <j!.2 = ±7r/2 

Note that the relative equilibria found here is identical to finding the equilibria of the second-order 
system given in Eqs. (17) and (18). The stationary conditions for energy also enforce that r = r = 4>2 = 
4>2 = 0- 

It is interesting to compare this result to the full case from Scheeres p3] when there is no assumption 
on the shape of the primary. In that case we have a dependence on <j)i , and the stationary conditions also 
tell us that (f>i = and 0i = 0, ±7r/2, or tt. Recalling the equation from Scheeres [19] . 
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The differences between Eq. (20) and (22 1 (aside from the normahzation) are that for the case when the 
primary is not exactly symmetric about its rotation axis, we must consider the entire angular momentum, 
Ktot, and therefore we also see Ji^ appearing in terms with J2- • Also, as the primary body becomes 
spherically symmetric, the dependence on the primary orientation disappears and — >■ Ii^ — Is. This 
causes the system to decrease the total number of equilibrium configurations from 8 to 4, and since all the 
dependence on the primary has been removed, the primary can also be spinning at an arbitrary speed, </>!, 
and the system can still be in a relative equilibrium. 

The location of the equilibrium points depends on three main parameters: the mass ratio, the angular 
momentum, and the body shapes. The influences of these three parameters are studied in Section [3. 5| 

Throughout the remainder of the paper, any reference to a nominal system refers to the system outlined 
in Table!!] which are based on the binary asteroid system 1999 KW4 [14, 7J. The new parameter that appears 
in Table 111 is 

X=^ (23) 

which is the ratio of the polar radius to the equatorial radius of the oblate primary body. This parameter 
becomes important for scaling the system. 



Table 1 The non-dimensionalized parameters for the nominal test system. 



V 


/3 7 


-Rp Re. X 






0.9257 


0.8109 0.6112 


2.3590 2.6182 0.9010 


2.4034 2.1175 


0.1973 0.2913 0.3434 



The nominal system has an angular momentum of K = 2.8382, which results in equilibria at r"*" = 9.2442 
with = —0.0497, and r" = 9.2869 with E~ = —0.0496 where the superscripts indicate which equilibrium 
point we are referring to; the + indicates the equilibrium point at 02 = 0, and the ^ refers to the equilibrium 
at <j>2 = 90°, as with . 



3.2 Equilibrium Point Linear Stability Analysis 



Given that the linearized system near the equilibrium point is simplified significantly, we can find the 
eigenvalues of the system analytically. The dynamics matrix around the equilibrium points has the form. 



A: 



10' 
1 
P Q 
i? S 



(24) 



where the non-zero partials (^, 
(P, Q, P, and 5"). The characteristic equation for this matrix is simply. 



and ) have been represented by simpler variable expressions 



(P + P + QS)y' + PR = 



The roots of this equation are. 



A 



2 _ P + R + QS± ^{P + R + QSy - APR 



(25) 



(26) 



The determination of the eigenvalues can easily be carried out numerically on a case-by-case basis to 
determine the dynamic behavior around the equilibrium point. 

The eigenvectors can also be computed analytically by writing the eigenvalue problem as. 



AI- A 







where the eigenvector is v = V2 V3 v^]^ . The four equations that result are, 

Xvi — V3 = 

XV2 — V4 = 



(27) 

(28) 
(29) 



6 



— Pvi + A«3 — QV4 = 

— RV2 = Svs + \V4 = 

The eigenvector can be computed from these relationships to be 

■ 1 ■ 
a 
A 
_a\ 

where 



QA 



In the nominal system, the eigenvalues of the equilibrium point r"*" are 



Al,2 

A3, 4 



± 0.0302i 
± 0.0362i 



(30) 
(31) 



(32) 



(33) 



(34) 
(35) 



which makes it a center point. However, we will refer to this equilibrium point as stable in the remainder of 
the study since trajectories can be bound around it, as will be explored. The eigenvalues of the equilibrium 
point at are 



Ai,2 = ±0.0307 
A3 4 = 0±0.0351i 



(36) 
(37) 



This equilibrium point has a unstable and stable asymptote associated with the real eigenvalues. It is 
interesting to note that the eigenstructure of this problem is nearly identical to that of the coUinear 
Lagrange points in the restricted three-body problem studied by Conley t5j. 



3.3 Equilibrium Point Energy 



It has been shown in Section [3. 1| that the equilibrium points are stationary points for the total energy. In 
this section, we determine if these stationary points are local maxima, minima, or saddle points for the 
energy. To find this, we must investigate the second derivatives of the energy with respect to the state, 
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where the subscripts indicate partial derivatives. 



Through investigation of Eqs. (123) - (125), it is quickly clear that, 

E, 
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The only cross second derivative remaining is with respect to r and (f>2. Looking back at the energy in 



Eq. (15), we see that the only part of the energy which contains both of these variables is the potential. 



Therefore, 



drd(j>2 drd(j>2 







(40) 



where the equality to zero at the equilibrium point was shown in Eq. (137). We have shown that all of 
the cross partials are zero, so that _Exx is a diagonal matrix. The diagonal entries of the matrix are the 
eigenvalues, and the definiteness of the matrix is found by looking at the sign of the eigenvalues. The 
diagonal entries are. 
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Err = V (43) 

E. . = (44) 

These results are straight forward to interpret. The velocity partials, Ef-f- and E^^^^ are always greater 
than zero. The angle partial, E^^^^, is positive for the equilibrium points at 02 = or tt, and is negative 
for the equilibrium points at <^2 = ±vr/2. The position partial, Err, is much more complicated and can be 
positive or negative depending on the system parameters and location of the equilibrium point. 

Combined, this tells us that an equilibrium point at (/!>2 = or tt can be a local minimum or a saddle 
point depending on the sign of Err- An equilibrium point at 02 = ±7r/2 is always a saddle point. The fact 
that the only energetically stable solutions are found at 02 = solutions was first shown by [l7]- In the 
nominal case, is energetically stable, and is a local minimum in energy. The other equilibrium point, 
r~ , is an energetic saddle, being a local minimum in the radial direction, but a local maximum in the 02 
direction. 



3.4 Osculating Orbit Elements at Equilibrium 



In considering a planar problem the orbital elements of interest are the semi-major axis, eccentricity, true 
anomaly, and argument of periapse. The behavior of these elements at equilibrium are discussed in this 
section. The other orbital elements, namely inclination and argument of the node are effectively meaningless, 
and won't be discussed here since they are either constant, undefined, and/or always zero. 
We can derive the Keplerian energy and angular momentum as. 
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where v is the inertial velocity vector. At equilibrium. 
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which was derived from Eqs. ( 12 ) and by setting ( 126 1 equal to zero as at equilibrium. The Keplerian energy 



and angular momentum at equilibrium are then determined as. 
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Then the osculating semi-major axis and eccentricity can be computed in terms of the energy and 
angular momentum as, 

-1 
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so that 

3 (li, -is + a 



2r2 



(52) 



Using Eqs. (50) and (51 1, we find that at equihbrium we have 

a = (53) 



Due to the relationship in Eq. (51 1, we can see that e > at the 4> = G equilibrium points because /i_ —Is > 
and C2 > 0. At the = 90 equilibrium points, e can be positive, negative, or zero because C2 can be 
negative. This tells us that at equilbrium the system is always locked at periapse or apoapse, depending on 
the sign of e. 



However, at a relative equilibria, the orbit rate is given by Eq. (471, and is constant. On a Keplerian 
orbit, this would imply that the orbit is circular, and this orbit rate is identical to the mean motion. In 
this case the eccentricity is in general non-zero and constant, and the semi-major axis is constant, along 
with the radius. Combining these results indicates that in fact the true/mean anomaly are constant (equal 
to or 180° as discussed above) and the argument of perigee is precessing at the orbit rate to enforce this 
condition. 

This can be shown by investigating the evolution of the eccentricity vector, which points to periapse. 
The eccentricity vector is defined by, 

e = V X - f (54) 



and in the secondary fixed frame at equilibrium the eccentricity vector becomes through use of Eqs. (47 1 



and (49) 



e = {r-^e'' - l)f = er (55) 



where e was defined in Eq. ( 52 I 



The rate of change of the eccentricity vector can be determined by using the transport theorem and 
resolving in the secondary fixed frame. At the </>2 = equilibrium point it becomes, 

e = eOy (56) 

and at the (f)2 = 90° equilibrium point it becomes, 

e = -edx (57) 

Given that in this rotating frame, the inertial rate of change of the unit vectors are 

i = 0y (58) 

ir = -e^ (59) 

we can see that the angle between an equilibrium point radius vector and the eccentricity vector is constant; 
therefore the anomalies (true, mean, and eccentric) are fixed at or 180°. If we assume that the argument 
of periapse is measured from a fixed direction in the x — y plane, then we can see that the rate of change 
of the argument of periapse is exactly equal to 9. 

In summary, a system which is at equilibrium will appear to an outside observer to be moving on a 
circular orbit. Due to the stability of the equilibrium points, this will likely only actually occur at the 
02 = point, which is commonly referred to as a synchronous orbit. If such a system is observed and fitted 
to Keplerian dynamics only, modeling each body as a point mass (or sphere) , the result would be a circular 
orbit with an eccentricity of zero. The computed semi-major axis would then imply an incorrect fi, and thus 
an incorrect mass of the system. Therefore we reiterate that it is crucial to account for the non-spherical 
shape of both bodies when fitting orbits to celestial objects. 

For reference and later comparison, the nominal system has a semi-major axis of a+ = 9.3269 and an 
eccentricity of e+ = 8.87 x 10^'^ at the <p2 = equilibrium point. The values at the (^2 = 90° equilibrium 
point are a" = 9.3270 and = 4.31 x 10"^. 
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3.5 Equilibirum Point Variation for Various System Parameters 



The basis for a system to have librational motion is in the properties of the equilibirum points. We have 
shown in the preceding sections that the equilibrium point at 02 = is spectrally stable and an energetic 
minimum; therefore this is the equilibrium point about with the system will librate. The equilibrium point 
at 4>2 = 90° is spectrally unstable and an energetic saddle, which means systems will generally not stay near 
this equilibrium point. Before moving on to study the actual trajectories of the system, we first study the 
locations and properties of the equilibrium points for different systems. The parameters that define these 
systems are the mass ratio, the angular momentum, and the shapes of the bodies. The effect of varying 
these parameters are discussed in this section. 

First, we investigate how the variation in mass fraction affects the equilibrium points. In order to isolate 
the effects from varying the mass fraction from the effects of the body shapes, we vary the sizes of the bodies 
along with the mass fraction to keep the moments of inertia constant as follows. Assuming equal density, 
the mass fraction is equivalent to the volume fraction, 

— ^1 _ R^Rp 

Vi+V2~ RlRj, + M ^ ' 

This can be solved for the oblate equatorial radius by using x to get. 

Using this equatorial radius, the moments of inertia of the oblate body become, 

\rI (62) 



5 

1 /„2 . „2\ 1 



Rl + Rl) = ^Rl (l + x') (63) 



' 5 

This paper is mainly concerned with systems where the ellipsoidal body is the smaller body, which 
would correspond to > 0.5. However, for the sake of completeness, we show in this section the locations of 
the equilibrium points for all values of v. When the ellipsoidal body is larger, our perspective changes and 
we think of this in terms of studying the location of orbits of an oblate satellite, instead of the libration of 
an ellipsoidal satellite. 

The locations of the equilibrium points are shown in Fig. [3] The general behavior is that as the mass 
fraction becomes smaller, the equilibrium points move to larger radii. The second plot shows the difference 
in radial distance between the two equilibrium points since they are too close to appear as separate lines in 
the scale of the first figure. The opposite trend is seen here in that the two equilibrium points are further 
apart radially as — >• 1. At the nominal value of angular momentum, the two equilibrium points are 0.04 
apart as was given at the end of Section [3. 1[ 

The astute reader will notice that this does not seem to match the results from Bellerose ^3j; this is 
because that work included a factor of v in their computation of angular momentum we do not include 
(e.g. Kg^iif,j.gi^f, = vK). The value of angular momentum used in this paper is absolute for any system. Also 
note that Bellerose discusses cases for which multiple equilibria appear; we find these cases as well for small 
values of K. However, since we are interested in studying the outer pair of equilibrium points that have the 
librational structure, we don't study the inner equilibrium points that may appear. 

The main result of varying the angular momentum as illustrated in Fig. [3] is that for higher values of 
K, the equilibrium points move to larger radii. This also has the effect of making the absolute difference in 
energy levels between the equilibrium points smaller, as is seen in Fig. |4] 

The shape of the bodies can have a number of effects on the equilibria. In terms of the oblateness of 
body 1, the more oblate the body, the larger the difference 7i, — Is becomes. This generally means that 
the location of the equilibrium points moves outward with a more oblate body. The effect on the 02 = 
and 02 = 90° equilibria are roughly the same. 

The effects of the secondary shape can be more varied due to the fact that the main asymmetry in this 
problem is due to the ellipsoid shape/moments of inertia. The effects are encompassed by the differences 
in and C2 , which can be changed due to variations of any of the three moments of inertia. Recall that 
changing also changes the value of the system moment of inertia, Iz- 

It is interesting to recall that the semi-major axis and eccentricity depend largely on the quantity 
(7i^ — Is + C^). While the equilibrium point at 02 = is always at periapsis, the 02 = 90° equilibrium 
can be at periapse, apoapse, or on a circular orbit depending completely on the moments of inertia of the 
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Fig. 3 The equilibrium location for varying values of the mass fraction are showrn. Three lines are plotted at different 
values of angular momentum: K = 2.8382, the nominal value in green, K = 3.9625 in red, and K = 2.1963 in blue. The 
dotted vertical line indicates the nominal value oi v = 0.9257. The lower figure shows the radial difference between the 
location of the two equilibrium points. 



bodies. Also, due to the same quantity, the 02 = 90° equilibrium point can become a local maximum in 
the radial direction (see Eq. (41)), although the point will still be an energetic saddle in general. 

In all cases studied here, the results of the stability of the nominal system equilibria holds for the outer 
set of equilibria points. Namely, the outer equilibrium point at <j!>2 = is spectrally stable and an energetic 
minima, while the outer equilibrium point at <j)2 = 90° is spectrally unstable and an energetic saddle. 
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2 3 4 5 6 7 8 

K 

Fig. 4 The diflference in energy between the stable and unstable equilibrium points, AE = Eun — E^tah^ is plotted in black 
for the nominal system. The three vertical lines indicate the three different values of angular momentum, classified by the 
same color scheme as used in Fig. [s] 
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4 Sufficient Condition for Bounded Motion 



When the system is not at equilibrium, the trajectory will evolve in r 
of motion given in Eqns. (17 1 - (181 



2 space according to the equations 
As this system is non-integrable, we are particularly interested in 



finding conditions which classify or restrict the trajectories that will occur. In this section, we present a 
sufficiency condition for bounded motion, which is taken to mean that the libration angle will be less than 
90° for all time. We then investigate how the sufficiency condition changes for different system parameters. 
Finally, we look at some properties of the system trajectories that meet the bounded sufficiency conditions. 



4.1 Bounded Motion 

The ZVCs can be used to investigate if the current system can become unbounded by reaching a phi = 90, 
which leads to the sufficiency condition for bounded motion: 

Theorem 4.1. Sufficient Condition for Bounded Motion 

Given a system with angular momentum K so that a stable equilibrium point exists at (j>2 = 0, the librational 
motion is bounded (4>2 < 90° always) if E < E^ . 



Proof: The relationship for the free energy of the system, Eq. (15), can be rearranged to become 



„2 12 



^ IK ^r, , \ 1 .2 \ lo vr 0).^ 



or equivalently, 



The right hand side of Eq. ( 64 1 is always positive, therefore we can state the relationship for a zero-velocity curve 
(ZVC), 

E-\^^^V{r,c^2)>Q (65) 

i?> ^;^ + V^(r-,02) (66) 

The stable equilibrium is the minimum energy location, and when E = E^ , the ZVC defines only the stable 
equilibrium. As E is increased, the ZVC will encompass larger areas mr — (f>2 space. If E is increased to E~ , the 
ZVC will touch the unstable equilibria, and because E~ is the minimum energy radius with <j)2 = 90° (see Section 



3.3), this is the minimum energy at which the ZVC allows the system to reach (t>2 = 90° 
That this is a sufficient, but not necessary, condition is shown in Section\5.1\ 



□ 

A given system will have some area in the r~(j)2 space which it can reside in based upon the free angular 
momentum and energy. Given a value for the free angular momentum of the system, the entire phase space 
can be mapped with varying energy levels determined by Eq. ( |66[ ). The easiest way to visualize this is by 
looking at the system from a secondary fixed frame, as shown in Fig. [2] For a given value of free energy for 
the system, there will be some area to which the primary is constrained to reside. Note that due to the fact 
that this relationship is an inequality, the primary can be anywhere inside the free energy level, not only 
on the surface. Therefore this relationship clearly doesn't solve the equations of motion to tell us what the 
state is at any given time, but it does tell us absolutely that the state is always inside the area bounded 
by that free energy. 

Theorem |4.1| is verified graphically for the nominal system in Figures [5] and [6] The nominal system 
zero- velocity curves are shown in Fig. [5] This is the typical structure for the zero- velocity curves seen in 
most situations where the libration between the two outer equilibrium points is being examined. In these 
types of plots, the stable minimum energy equilibrium point is along the x-axis {(1)2 = 0), while the unstable 
equilibrium point is on the y-axis (02 = 90°). As the energy is increased from the minimum at the x-axis 
equilibrium point, the zero- velocity curves allow for larger libration angles until the unstable equilibrium 
energy is reached. At energies above the unstable equilibrium point energy the ZVCs open (becoming two 
separate ZVCs) and the secondary is free to circulate in this area. The black shape around the origin is 
the projection of the ellipsoid plus the equatorial radius of the oblate body; if the energy is high enough so 
that this region is inside the ZVC bounds then an impact between the bodies is possible. The color bar lists 



the values of energy {E) corresponding to each ZVC. Recall from Section 3.1 that the stable equilibrium 
point has energy E^ = —0.0497 and the unstable equilibrium has an energy of E^ = —0.496, where the 
difference between the two is 5E = 1.500 x lO""'. 
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Fig. 5 Phase space for the nominal test system described in Table jl] The closed ZVCs are at increments of energy of 10% 
of SE, so that the first close curve has energy E = _E+ + 0.1<5i?. This first curve limits the libration angle to a maximum 
of 02 = 18.4°. 



In order to make clear the behavior of the energy in the vicinity of the equilibrium points, we plot 
what is effectively a cross section of Fig. [5] in Fig. [6] This clearly shows the variation of the energy in the 
radial and circumferential directions. It is clear that both of these equilibrium points are minima in the r 
direction, however only the x-axis equilibrium point is also a minima in the (j>2 direction. 

It should be noted that the same behavior is seen around the stable equilibrium point at 02 = 180°, and 
Theorem |4.1| can be applied there as well. In this paper, we generally only look at the 4>2 = equilibrium 
point for clarity. 

Consider a zero velocity curve with the free energy, E, and the free angular momentum, K . When the 
inequality in Eq. (65) is precisely equal to zero, we know that f = and 02 = 0; all of the free kinetic 



energy in the system has been transferred to the potential energy. At any point inside the zero velocity 
curve(s) defined by E, the inequality will be greater than zero as there is an excess of energy defined from 
the ZVC definition, 



AE{r,(t>2) = E- 



2 h 



V{r,cj>2) = -mr 



(67) 



This situation is depicted in the cartoon shown in Figure [Tj Returning to Eq. ( 64 1 , it is clear that in this 
case, the left hand side is greater than zero, and therefore either jfj > 0, \<i>2\ > 0, or some combination of 
the two. 

Note that the excess energy, AE is not constant on a given trajectory as it varies with the position 
variables. The maximum value that AE can reach is at the 4>2 = equilibrium point since this is the 
minimum potential energy location. At any given point on a trajectory, the excess energy can be used to 
determine the range of possible values of the velocities based on the kinetic energy in the radial and spin 
components. 
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Fig. 6 Variation of the energy near the equilibrium points for the nominal system. The upper plot shows the variation 
with respect to the radial distance at <f)2 = and <j)2 = 90°. The location of the equilibria are marked by the vertical dotted 
lines. The lower plot shows the variation in the energy at the stable equilibrium point radius as (f>2 is varied. The points 
at ±n/2 on this plot are actually equivalent to where the blue vertical dotted intersects the red line in the upper plot, not 
precisely at the unstable equilibrium radius. However, the same trend of being a local maximum with respect to variation 
in 02 holds at the unstable equilibrium point radius. 



The split in the kinetic energy wiU be determined by the factor k, which dihneates what percentage of 
the excess energy goes into the radial velocity. The scale factor is therefore bounded such that 

< K < 1 (68) 

This means that the radial velocity and ellipsoid rotation rate are defined as, 

= ^ (69) 
u 

= ^(i^^ (70) 

Fig. |8] shows the trajectories of three cases for the nominal system with a closed ZVC, each beginning 
with a different amount of excess energy in the radial and angular forms of the kinetic energy. These results 
demonstrate that the ZVC is indeed accurate, and that the trajectories stay bound within this area of 
r-4>2 space. Furthermore, depending on the initial velocities, the area within the ZVC that is filled by each 
trajectory is limited further; the trajectory with all the initial velocity in the radial direction (plotted in 
blue) explores the phase space the furthest in the radial direction, but not as far in the angular direction. 
The opposite is true for the case with all of the initial velocity in the angular direction, which is plotted in 
red. 

The ZVCs are a crucial property of the system which limit the phase space in which the system can 
reside. Since the conditions in Theorem |4 . 1 1 are so important, we investigate how the ZVCs vary when the 
system parameters are changed. 
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Fig. 7 The zero-velocity curve for the given energy, E, has zero excess energy, as shown. Moving in toward the equihbrium 
point increases AE until reaching the equihbrium point (marked by the x), where the excess energy is at a maximum. 
Intervening zero- velocity curves are drawn with dotted lines. 
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r 



Fig. 8 Three trajectories for the nominal system with an energy oi E = + 0.3SE, which means the zero-velocity curve 
is closed, as shown in black. All three trajectories are initiated at r"*" with k = 0, k = 0.5, and k = 1. 
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4.2 Parameter Influence on Zero- Velocity Curves 



Eq. ( 65 1 shows us that the zero- velocity curves are defined first and foremost by the angular momentum 
and energy in the system. However, the exact form of the zero-velocity curves depends on all the system 
parameters such as the moments of inertia of the bodies, and the mass fraction of the system. In this section 
we explore how the nominal zero- velocity curves shown in Fig. [s] change when the system parameters are 
varied. 

4-2.1 Variation Due to Mass Ratio 

In this section, we look at the ZVCs for three different mass ratio systems pulled from Section |3.5[ By 
comparing these different mass ratio systems to our nominal system in Fig. |5j we get a good idea of the 
effect of changing the mass ratio on the ZVC structure. 

The first case is the nominal system scaled to a mass ratio oi v = 0.5, shown in Fig. [9] Even with the 
scaling of the bodies to equal masses, the structure of the ZVCs is very similar. The values of the energy 
and the location of the equilibrium points have changed due to a change in angular momentum, but the 
interesting point is that the behavior of the system between and around the equilibrium points found in 
our nominal case of a relatively small ellipsoid orbiting a larger oblate body appears to hold all the way to 
equal masses. 




1 2 3 4 5 



Fig. 9 Zero- velocity curves for the scaled system with v = 0.5 and K = 1.0007. The first 10 zero-velocity curves surrounding 
the stable equilibrium point are plotted in increments of energy of O.ISE. 



The second case we look at is with a mass ratio of v = 0.05 shown in Fig. |10[ This system was chosen 
because it illustrates a case discussed by Bellerose 0] where there are two equilibrium points on the x-axis, 
and the outer point is spectrally stable and an energetic minima, while the inner point is spectrally unstable 
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and an energetic saddle. It turns out for these systems that there are also two equilibria on the y-axis. As 
expected from the systems studied so far, the outer y-axis equilibria is spectrally unstable and an energetic 
saddle. It is interesting to find that the inner y-axis equilibria is spectrally stable while being an energetic 
maxima. 

This is interesting, as it implies that as v ^ 0, so that the system being studied is that of a small 
sphere/oblate body orbiting a large ellipsoid, there will be stable orbits along the y-axis around three times 
the ellipsoid semi-major axis. The orbits that have the small body near the same radius on the x-axis will 
be unstable. 




2 4 6 8 10 

X 

Fig. 10 The zero-velocity curves for the system with u = 0.05 and K = 0.1479. The minimum energy point is the outer 
equilibrium point on the x-axis, however low energy trajectories can exist inside the inner pair of equilibria, although they 
would likely quickly impact. 



The final case we show is that with a mass ratio oi u = 0.15, shown in Fig. |11[ The y-axis equilibrium 
points keep the same general behavior as those in Fig. |10[ however in this case the x-axis equilibria have 
vanished. This appears to make low energy impacts possible between the two bodies starting from basically 
anywhere in the system; in other words the energy decreases all the way to the surface along the x-axis, 
and since a given trajectory can explore the region between ZVCs for its energy, they can reach the surface 
along the x-axis. 

In the cases in Figs. [lO] and |11[ it is important to note that there may be some non-trivial errors in 
the values obtained for the closer equilibria due to the fact that we are using a second order expansion 
for the potential of the bodies. If accurate results are needed for these cases, we suggest calculating the 
outcomes with the exact potential expression used by Bellerose [3]. However, the qualitative results match 
these previous investigations with the exact potential. Furthermore, with the main interest of this paper 
being the outer set of equilibrium points, these points are far enough from the bodies that the differences 
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Fig. 11 The zero-velocity curves for the system with v = 0.15 and K = 0.3211. At this configuration, the x-axis equilibria 
have disappeared and the path to impact along the x-axis is open for all energies. 



between our second order approximation and the exact potential should be small enough to make little 
difference. 
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^.2,.2 Variation Due to Angular Momentum 



As was alluded to in Section [3. 5 [ the main effect of changing the angular momentum is to move the equilibria 
in or out. This effect is also seen on the ZVC structure in Fig. [12] the ZVCs have the same structure as the 
nominal case, they are basically shifted in or out. It was demonstrated in Section [4. 2. 1| that a combination 
of low 1/ and K can cause a second set of equilibrium points to appear, however this situation doesn't 
happen for larger values of u due to the fact that these solutions to Eq. ( 20 1 are less than the combined 
radii of the bodies, and usually less than the radius of the secondary. 




Fig. 12 The ZVC structure for two cases of the nominal system with different angular momentum values. The case on the 
left has a lower angular momentum than the nominal a,t K = 2.1963, while the right figure has a higher value of K = 3.9625. 
Recall that the nominal angular momentum is K = 2.8382. The main effect of changing the angular momentum is to move 
the equilibria out for higher values, and in for lower values. 
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4-2.3 Variation Due to Body Shapes 



The degree to which the secondary is triaxial plays a large role in the determination of the zero-velocity 
curves. This is investigated by tweaking the nominal case with two different values of the intermediate 
moment of inertia, 12^, first so that it is just greater than the /2^, and second so that it is just less than 
■ The main contribution of modifying the ellipsoid moments of inertia, however, is that the difference 
{l2y — l2^) controls the dependence of the potential and therefore the ZVC on the libration angle. The 
ZVCs for these two cases are compared to the nominal case in Fig. [13] The ZVCs of the same color are at 
the same energy across the three figures, making it easy to see that as the intermediate moment of inertia 
is made smaller, the body can circulate much easier. 




Fig. 13 Zero-velocity curves for the nominal system with three different values of the ellipsoid middle moment of inertia. 
From left to right, the minimum value is l2y =, the nominal value is l2y =, and the maximum value is l2y =■ Each of the 
three figures has ten ZVCs plotted at the same energy levels, clearly illustrating that it is easier to cause an ellipsoid with 
a smaller difference {l2y — 1 2^), the minimum case, to begin circulating. 



Note that, as discussed in Section [3. 5[ changing the oblateness of the other body will move the location 
of the equilibrium points in or out. However, there is very little difference beyond this radial shift on the 
ZVCs, so this effect isn't illustrated. 
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4.3 Time- Variation of Osculating Orbital Elements 



The orbital elements of a non-equilibrium trajectory are no longer useful as metrics for the system because 
they change constantly throughout the course of the trajectory. The groundwork for calculating the orbital 
elements has already been laid in Sections |3.4| and |4.1[ At any given point in r — 02 space, given the energy 
and angular momentum, the possible values of r and <^2 can be calculated using Eqs. (701 and (691. The 



Keplerian energy and angular momentum can then be calculated from Eqs. (461 and (451, which in turn 
allows us to calculate the semi-major axis and eccentricity with Eqs. (501 and (51 1. Note that for these 
calculations, the sign of r makes no difference, however the sign of 02 directly effects the value of 9, which 
then changes the value of the Keplerian energy and angular momentum. 



The eccentricity vector can again be calculated through Eq. (54). Note that in this case, r is no longer 
zero in general, so that 



(71) 



where 9 = x f. This implies that whenever f ^ 0, the eccentricity vector no longer lies along the f 
direction, and therefore the system is no longer at an apse. 

To illustrate the behavior the semi-major axis and eccentricity on a non-equilibrium trajectory, the 
corresponding values of the semi-major axis and eccentricity for each of the trajectories in Fig. |8] are shown 
in Fig. |14[ As with the radial-angular phase space, each trajectory occupies a distinct region of a — e space 
despite the fact that they each have the same energy. 
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Fig. 14 Eccentricity versus semi-major axis for the three trajectories shown in Fig. [s] 



To further understand the possible values of semi-major axis and eccentricity that can be reached along 
the trajectories, the values in the entire region inside a give ZVC can be mapped. Figs. |15| - 116| show these 
results. Each subplot can be imagined as a cross section of a ZVC plot (e.g. Fig. ItI) at a different value 
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of cj>2, where the out of plane axis teUs us how the excess energy is apportioned between the radial and 
angular velocities through the use of the parameter t], which is defined as 



f - K 02 > 



1 



< 



(72) 



so that negative values of 77 correspond to negative values of 02 • The radial bounds of each subplot corre- 
spond to the inner and outer edges of the ZVC for this energy level, which is why at larger values of 02, 
the radial range is smaller. 

From these figures, we see that the semi-major axis reaches its maximum and minimum values near 
the stable equilibrium point, with all of the energy in negative or positive angular velocity, respectively. 
The eccentricity, on the other hand, reaches its maximum value near the inner boundary of the ZVC at 
02 = with most of the energy in the radial velocity, but some in a negative angular velocity. The minimum 
eccentricity, which is zero, is seen just outside of the equilibrium point distance with most of the energy in 
either positive or negative angular velocity. Note that the colorbars near the last subplot are valid for all 
values of 02. Clearly a wide variety of combinations of semi-major axis and eccentricity can be seen over 
the course of trajectory. 





Fig. 15 Semi-major axis values for various libration angles with the same energy as 



Fig.js] 
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5 Transient Dynamics 



In this section, we explore trajectories for systems with open ZVCs (those with E > E^). Specifically, we 
explore the existence of trajectories that behave as if they are bounded, but do not meet the sufficiency 
conditions outlined in Section |4j These results prove that Theorem |4.1| is in fact a sufficient, but not 
necessary, condition for bounded motion. We also investigate the special case of bounded motion for periodic 
trajectories, and the families of periodic trajectories that exist as the energy is varied. 



5.1 Existence of Bounded Motion with Open Zero- Velocity Curves 

The existence of bounded trajectories with open ZVCs is shown through numerical examples. Two such 
trajectories are shown in Fig. [17] Although the ZVC is open, not every trajectory circulates. This is an 
extension of what was seen above with the closed ZVC trajectories; depending on the initial conditions of 
the trajectory at (l>2 = 0, the trajectory will only explore a portion of the r — (j>2 phase space inside the 
ZVC. The existence of this bounded trajectory proves that Theorem |4.1| is a sufficient, but not necessary 
condition, as this example does not meet the requirement of the theorem. Furthermore, using the ZVC 
corresponding to the energy of the unstable equilibrium point is a conservative estimate of the energy at 
which a system will have a circulating, as opposed to librating, secondary body. 




5.2 Periodic Orbits 

Periodic orbits can be found in this system when a given orbit repeats itself exactly after a certain period. 
These orbits are found when the state transition matrix or monodromy matrix have unity eigenvalues. A 
detailed methodology for finding periodic orbits in this manner is given in I2- In this system, due to the 
symmetry about 02 = 0, most of the periodic orbits will have a crossing of this line of symmetry with 
f = 0, although this is not strictly required for a periodic orbit. However, it is true that every trajectory 
that crosses 02 = twice with f = is periodic. The period may be very long, but due to the symmetry of 
the problem it will eventually repeat. 

Therefore, our search consisted of starting at different radii with f = 0, and looking at subsequent 
crossings of 4>2 = 0. The results of this search are combined in a Poincare map at of the r ~ r space 
at 02 = 0. One illustrative case is shown in Fig. [18] where the crossings of several periodic orbits are 
highlighted. 

The actual trajectories of several periodic orbits in the r— 02 space are shown in Fig.|19[ corresponding to 
the highlighted cases from Fig. [18] As alluded to previously, there are actually an infinite number of periodic 
trajectories in this system, with most of them having very large period ratios such as the T = 49.04 and 
T = 56.84 orbits shown. Of particular interest are the two periodic orbits with near T = 1. In the example 
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Fig. 18 Poincare plot of r and r for the nominal system with E = i?+ + 0.3SE, the same as in Fig. [s] Five different 
periodic orbits are highlighted. The small blue dots correspond to many other integrated trajectories that are not periodic, 
and so they pierce the Poincare map at different locations each time. 



shown in Fig. |19[ these are the T = 0.9896 case, which orbits counter-clockwise, and the T = 1.274 case 
which orbits in a clockwise direction. These two orbits roughly span the space inside the closed ZVC in 
this they are associated with the eigenmodes of the equilibrium point. 

The evolution of the eigenmode periodic orbits with changing values of energy are explored in Figs. [20] 
and [21] The first plot shows many different trajectories as the energy is increased. Particularly interesting 
behavior is seen with the clockwise trajectories. When the energy is approximately 75% of the way to the 
unstable equilibrium point energy, the shape of these orbits changes from the oval seen in Fig. [19] to a bow 
shape, and eventually (at the highest energy levels shown) the trajectory crosses itself near 02 = ±50°. 
The maximum libration angles reached by these periodic orbits is also very high, approaching 80° at the 
energy levels shown. By comparison, the evolution of the counter-clockwise orbits is much more regular. 
The shape stays basically the same throughout, simply growing in size as the energy grows. Note that 
the maximum amplitudes reached in these cases is much smaller than the clockwise orbits, getting to a 
maximum of roughly 02 = ±25°. 

A different view of the trajectories is shown in Fig. |21[ This plot shows the radial locations of the 
02 = crossings with f = at each energy level. This clearly shows the more regular evolution of the 
counter-clockwise orbits. It is especially interesting that these periodic orbits exist well beyond the energy 
of the unstable equilibrium point, as this gives clear examples of stable, libration bounded non-equilibrium 
orbits with open ZVCs. 



27 



T = 56.84 




r 



Fig. 19 The trajectories for five distinct periodic orbits in tlie r-02 phase plane which were highlighted in Fig. |l8| The 
trajectories are labeled by the ratio of their period with the period of the equilibrium orbit. Note that the T = 56.84 
trajectory was plotted in cyan in Fig. |18| and the T = 49.04 was plotted in magenta. 




Fig. 20 The periodic orbit trajectories a re s hown, with the clockwise orbits (T = 1.274 in Fig. |l9| in the left plot and 
the counter-clockwise (T = 0.9896 in Fig. \l9\ orbits in the right plot. For this test, the energy was varied from E = i?+ 
to E = E~ + 5E. Note that the lower energy cases are shown in dark blue, varying through green and yellow, up to the 
highest energy cases in dark red. 
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Fig. 21 Radial locations of tiie (f>2 = crossings with r = at varying energies. The clockwise orbits are shown in red, and 
the counter-clockwise are shown in blue. For reference, the stable equilibrium point energy is shown as a vertical dashed 
line, where the periodic orbit families disappear. The unstable equilibrium point energy is also shown as a dash-dot vertical 
line. 
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6 Sufficient Conditions for Unbounded Motion 

In Section [4] we presented a sufficient condition for bounded librational motion. In Section [5j we proved 
this condition was only sufficient due to the existence of bounded trajectories when the conditions from 
Theorem |4. 1 1 are not met. Significantly there are families of periodic orbits that exist even when the ZVCs 
of the system are open. In this section, we analyze the converse problem of determining when unbounded 
librational motion occurs. When the system has a mass fraction near unity, conditions for unbounded motion 
can be derived analytically. For non-unity mass fractions, strict conditions are unavailable. However, we 
present analysis for these systems that illustrates the structure of the problem and initial conditions that 
will lead to circulation. 



6.1 Analytical Limits for Systems with u ~ 1 

In this section, we show sufficient conditions for circulation for systems with mass ratios close to unity. 
This restriction limits the effect of the coupling on the orbit of the small ellipsoidal body, which allows 
us to assume that the orbit itself remains unchanged due to the libration of the secondary. However, due 
to this assumption the results can only be claimed as sufficient conditions for circulation, as shown in the 
following theorem. 

Theorem 6.1. Sufficient Condition for Circulation 

Given a system with i/ ~ 1, so that the orbit can be considered unperturbed by the ellipsoidal secondary, with some 
initial secondary spin rate (j>2.o oi a libration angle of (f/ = 0° at some location on the orbit (ro), the secondary 
body IS guaranteed to circulate if the spin rate satisfies the relationship 



12 ^ h,0 
»2,0 > = 
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where ra is the apoapse radius of the orbit. 

Proof: The free energy of the system is written using Eqs. Q, ([9| and (10) as, 
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where -Eo^fc orbit energy that is independent of the librational state, E^ is the librational energy (identical 

to a pendulum), and Ecoup is the "coupling" energy. The first two energy terms are given by relationships 
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In the case of a point mass, the orbit energy collapses to the Keplerian energy of the system. Recalling the 
trigonometric relationship 
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allows us to determine the coupling energy from Eq. (74 1, using Eqs. (151, (751 and (761 as 



(77) 
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where C2 was defined in Eq. ( 2 1 1 . 

// the orbit energy is constant, then the available energy is determined as 

^Eiii, = E — Eorb = E4, + Ecoup 

Since this value is a constant at any point on the trajectory, we can relate two points by 

^Eiihiro, (1)2 = 0, 02,o) = ^Eiibirm, (t>2,m, <i>2 = 0) 



(79) 



(80) 



This expression relates the available energy for a point on the orbit at radius rp with some initial libration rate 
02.0 and zero libration angle, to another point on the orbit at radius rm with libration angle 4>2,max and a libration 
rate of zero. In other words, this tells us the maximum achievable libration angle at radius rm given a system 
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that had an initial lihration rate at rg. Note that this relationship is independent of r. Writing out the available 



energy at these two locations using Eqs. ( 76 I and ( 78 I gives 

ri2 



21. 



z,0 



3u 

m. 



(■f2„ - 1 2,) Sin^ </>2,m - ^ (/l. 



To determine the bounding condition, we look for the case when 02,™ = 90°, and solve for (j>2,o to get 

Iz,0 



'i'2,0 



> 



I2^+IU 



1 



(81) 



(82) 



where the fact that 



-2/2, + /2„ + I2, was used. 



For any given starting location, the highest 02, viill be required if rm is as large as possible, so rm = ra- 
Therefore if the the libration rate is high enough to circulate the secondary at apoapse, then the body will circulate 
no matter what point on the orbit it reaches its maximum libration angle. Thus the bounding condition to ensure 
circulation occurs when 



;2 ^ Iz,o 
92,0 > = ^ 



+ (12 



+ I1 



r-3 



Theorem |6 . 1 1 assumes that the orbit perturbations caused by the hbration of the secondary are ignorable. 
Therefore the orbit can be described using Keplerian orbital elements, which will be constant. We can 
express the initial orbit radius in terms of the eccentricity, semi-major axis, and true anomaly as 



and then the apoapsis radius is 
and finally the periapse radius is 



1 + e cos / 
a(l + e) 



-0 = (83) 



(84) 



rp = a{l - e) (85) 
Using these relationships, the condition for circulation can be expressed in terms of the orbital elements 



"2,0 



> 



1 + e cos / 
a(l-e2) 



+ {I2, - l2^ + lu 



1 - 



[l^ef 



{l + ecosfy 



(86) 



Theorem |6.1| gives the sufficient condition for circulation for a given starting point, however we can 
define a single condition for circulation that is sufficient for any location on the orbit by defining the 
highest value that leads to circulation, which is referred to as the uniform condition. 

Remark 6.1. Uniform Condition for Guaranteed Circulation 



f 



The highest bound for circulation is given by maximizing Eg. (86) in terms of f. The maximum occurs when 
0, or ro = 



and the circulation condition becomes 
I. 



"2,0 



> 



/2.a5(i_e)5 



(/2„-/2j + (/2. -/2. -/.) 



2e (3 + 



(87) 



This condition is considered uniform because it combines the extremes for both ro and rm in terms of 



maximizing the required 02,0- The condition given Eq. (87) is the baseline sufficiency test for circulation. 
At any point on an orbit, if 02,0 meets this condition with 02 = 0°, it is guaranteed to circulate directly, 
meaning that the current trajectory will pass through ±90° before returning to 02 = 0°. 

Depending on the actual location on the orbit, the value of 02. which is sufficient to cause circulation 
will generally be lower, and is given by Eq. ( 73 ). It is important to note here that Theorem 6.1 is a sufficiency 



condition only, but it is possible that a lower value of 02, will lead to circulation. This is because Theorem 
6.1| depended on the most extreme value for Vni- However, due to the fact that it is easier to circulate at 
all rm < ra, it is true that lower libration rates will lead to circulation in many cases. This depends on the 
phasing of the librational motion with the orbit to determine at what radius cj>2,max is achieved, and if 02, 
was large enough at the 02 = crossing preceeding this 02,ma2;, 

circulation can occur. Therefore the lowest 
possible spin rate that can lead to circulation can be determined in the opposite limiting case where 02, 
is minimized with respect to rm and ro. In fact, this becomes a sufficient condition for bounded libration. 
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Corollary 6.1. Sufficient Condition for B ounded Motion with ~ 1 

Given the system considered in Theorem 6.1 the lowest value of (j>2.o which can lead to circulation is computed 



so that any value below this will not lead to circulation 



„3\ 1 

^0 



where rg is found as the real positive root of, 
'■o 



31^4 ( -Is+ C+ 



2/2 



2 
^0 



2(/2 



+ h 



(89) 



if the resulting < ra, otherwise rp = ra- 

Proof: Eq. \8S\ is minimized with respect to ro and rm ■ The partial with respect to rm shows that the function 
monotonically increases with rm, therefore the minima with respect to rm «s at the lower boundary, rp. 

There exists an extrema of Eq. 82 with respect to ro, which is located at the root of Eq. (89). Using the Routh 



criteria it can be proven that there is one and only one real positive root. It can also be shown that this extrema 
is in fact a minima by investigating the second partial with respect to rg . However, the location of this minima 
depends on the system parameters, and it is not guaranteed to lie within the range of possible radius values given 
by rp and ra. Therefore, if the minima is found to be greater than ra, the constrained minima is at rg = ra. 



These results outline a set of bounds that can be used to determine if a system will circulate or not 
depending on the libration rate at 02 =0°. The main boundaries to be tested to give information for an 
entire orbit are given by Remark |6.1| and Corollary |6.1| If either of these conditions are satisfied, we can 
immediately state if the body will or will not circulate. 

When the libration rate is between these two boundaries, further inspection must be made. Theorem 



6.1 can be used to determine if a given initial condition (rp and 02, o) will circulate directly. However, this 
relationship does not tell us if the body will ever circulate; only if it will circulate directly. The difficulty 
is in the phasing of the libration and the orbit. If the initial libration rate is below the bound given in 
Theorem |6.1[ then the body will not circulate this on this oscillation. However in general the body will pass 
back through 02 = 0° at a different state, and then the condition must be checked again. 

Corollary |6.1| is an extension of Theore m |4.1| for the specific case when ~ 1. If the total energy in 
the system is low enough so that Theorem |4.1| is applicable. Corollary |6.1| will also be satisfied. However, 
Corollary |6.1| can be satisfied when Theorem |4.1| is not applicable. 

If circulation does not occur, it is useful to be able to determine the libration amplitudes that can be 
expected. This is addressed in the following corollary. 



Corollary 6.2. Maximum Libration Angle 

Given the system considered in Theorem 
angle that can be obtained is given by 



6.1 



1/02,0 is below the limit for circulation, then the maximum libration 



1 



3(^2„-/2j 



+ C+ 



(90) 



Proof: This condition is determined by rearranging Eq. (81) to get 



sm 



P2,max 



3(^2, 



'2.0 



+ h 



'z,0 



+ C+ 



1 - 



(91) 



This relationship is maximized with a minimum rm: giving rm — the periapse radius. 



If the maximum libration relationship from Corollary 



is tested with 02, o above the circulation limit, 
the right hand side of Eq. (901 will be greater than one, resulting in an imaginary value for 4>2,max. This 



indicates the body will circulate. The maximum libration angle that would be reached at any point on the 



orbit can be determined by using a radius other than the periapse radius in Eq. (91) for rm 
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Remark 6.2. Simplification to Classical Results 

Classical gravity gradient results assume that the orbit is circular (rm = = r) and that the system is domina ted 
by the orbit so that ur^ ^ l2^, and therefore Iz — >■ vr"^ , and u = 1. Substituting these conditions into Eq 
gives the classical result 



81 



2 3 M2„ - ^2^ \ . 2 



l>2 



r3 



which is recognized as the energy integral of the pendulum equation of motion 



sin (j}2,max (92) 



3 / -^2„ - 

2^ 



sin 2(^2 (93) 



Derivation of the classical results are given in [2] and 9_ among many other sources. Note that Eq. 
( 93 1 has a slight difference to the referenced results because we have normalized the units so that /i has 



disappeared. 

These results can be illustrated through several numerical systems. First, a basic example is shown in 
Fig. 22 with very low eccentricity (e = 2.4 x 10~®). The full equations of motion given in Eqs. (17 1 - (181 
are integrated for two cases, where one trajectory is initialized with 95% of the excess energy is given to 
4>2, while the other has 70% apportioned to the initial libration rate. Both cases start at r = r+ and f = 0. 
Using Theorem |6.1[ we can see that the former case is above the limit and therefore will circulate, while 
the latter case is below. The maximum libration angle is calculated for the k = 0.3 case from Eq. (91 1 with 
r-m = r'^ , and is found to closely agree with the maximum libration amplitudes seen during the simulation. 



zvc 




5r 

Fig. 22 Two different trajectories for tlie nominal case with energ y 10 % above tlie unstable equilibrium energy. The dotted 
horizontal line indicates the computed <j)2,max = 61.34° from Eq. ( |90| l for the k = 0.3 trajectory. Note that on the x-axis, 
<5r = r — r"*". 



The results developed in this section can also show when the excitation of the libration amplitude from 
the eccentricity of the orbit can cause the secondary to begin circulating. Application of Corollary |6. 2| for a 



range of ro values is shown in Fig. 23 In this case, every case is tested with 02,0 = 02, o = 0, so the excitation 
occurs only from the presence of available energy for this orbit due to variation in the radial position. The 
results are presented in terms of eccentricity and true anomaly, which determine the extents of ro, but are 
more intuitive for this application. It is shown that for this case, if the eccentricity is less than 1 x 10~^ 
the libration is bounded for any initial condition. For larger eccentricities, the initial conditions must be 
nearer to periapsis in order for the libration to remain bounded. It should be noted that each eccentricity 
line is at a different energy level, but they all have the same angular momentum value for this study. 
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In order to illustrate the validity of these relationships, four different simulations indicated on Fig. [23] 
with circles were simulated for 500 orbits. These trajectories are plotted in Fig. |24[ It is clear that the 
bounds are accurately predicting the librational motion, and the unbounded case does circulate directly as 
predicted by Theorem |6.1| for that case. 
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Fig. 23 Maximum libration amplitudes at different points on a Keplerian orbit with ip2 = 0. When e = 0, maximum 
libration amphtude at apoapsis is only 0.34° and at e = 1 X 10~®, the amplitude can reach 1.65° at apoapsis. The unstable 
eccentricity is between e = 1 X 10"** and e = 2 X 10~*. For all eccentricities, there are some initial conditions that stay 
bounded; e.g. for e = 0.05, the libration amplitude is bounded for u < 0.0125°. 
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Fig. 24 Four trajectories chosen from Fig. |23[ integrated for 500 orbits. Integration of the circulating trajectory was 
stopped upon passing </<2 = 95°. These numerical integrations verify the validity of the relationships shown in Fig. |23| 



It is interesting to investigate the case of Saturn's moon Hyperion, which is known to be in chaotic 
rotation. Using the moments of inertia for Hyperion from ^ (renormalized by a = 180 km for Hyperion), 
and the known Hyperion orbit and Saturn properties, Eq. (87) is used to find the <f)2 that will guarantee 



circulation. We find that if Hyperion were to have |<^2l > 0.00524 degrees per second at periapse with zero 
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libration angle, it will enter a circulating state. Given the extremely small limit found, it is no surprise that 
Hyperion is in an unbound chaotic rotation state. 
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6.2 Initial Condition Mapping to Circulation for the General Case 

In the general case with u < 1, the relationships developed in Section [6. 1| can't be used because the orbit 
energy can not be considered constant due to the coupling from the secondary body's motion. In this 
section, we investigate which initial conditions at cj>2 = circulate numerically. Our approach is to sample 
the possible initial conditions at 4>2 = 90° and propagate them backwards in time to see where they originate 
at 4>2 = 0. 

By fixing one of the states (we choose to fix 02), the remaining states are required to lie on the surface 
of an ellipsoid of constant energy. We use these constant energy surfaces to investigate the states of the 
simulations at 4>2 = and 4>2 = 90°. These visualizations are essentially 3-D Poincare surfaces. The case 
investigated here is for the nominal system with E = + 0.15 E. 

The constant energy surface at 02 = 90° is illustrated in Fig. [25] Initial conditions which sampled the 
entire surface were tested, with those that propagated back to 02 = shown as black dots on the ellipsoid. 
There were 3362 initial conditions tested, exactly half propagated to cj>2 = 0, the other half propagated 
toward 02 = 180° and were ignored for this analysis. As seen in Fig. |25[ the split is basically associated 
with the sign of 02; those with a positive libration rate when crossing 02 = 90° generally came from 02 = 0. 
However, it is interesting to note that there 22 of the 1681 cases did have negative libration rates, seen as 
those will small values of r on the ellipsoid. 




r 

Fig. 25 Constant energy surface at 02 = 90° for the nominal case with E = E~ + 0.1<5_E. The black dots indicate the ICs 
tested which, when propagated backwards in time, pass through 02 = 0. 

The evolution of the backwards propagated trajectories are shown as they cross 02 = in Fig. [26] The 
initial point where the initial conditions reach 02 = are labeled as the zeroth crossing, and plotted also as 
black dots. We then plot each of the next six crossings in different colors, alternating between the north and 
south poles of the ellipsoid as the trajectories cross in opposite directions. Note there is definite structure 
to the region of each crossing on the ellipsoid. 

Fig. |27| shows the same thing as in Fig. |26[ except that all crossings after the sixth are plotted with red 
dots. This shows that although the coverage continues to spread and lose some of the distinct structure of 
the earlier crossings, there are still interesting features. First, there is a large area surrounding the equator 
where there are no trajectories that are connected to circulation. This tells us that trajectories that reside 
in this region will stay bounded. We point out that the points around the north pole are biased as a group 
toward the smaller radius values, while the south pole cases are biased toward the larger radius values. 
Second, there are some significant areas on the ellipsoids amongst the propagated trajectories that are not 
filled, such as on the north pole view around f = and small values of r. These gaps correspond to the 
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-0 02 8-6 r f 10 02 

Fig. 26 The north pole (left, positive values of <j)2) and south pole ( righ t) of the constant energy surface at (/>2 = for the 
nominal case with E = E~ + 0.1<5-E. The initial conditions from Fig. |25| are integrated back in time, with each of the first 
seven crossings of the (/)2 = surface plotted in different colors as indicated. 



periodic trajectories discussed in Section |5.2| The third note about this plot is that currently the region 
immediately surrounding the south pole is unpopulated. These states would correspond to trajectories that 
would evolve through 02 = —90°, as opposed to through 4>2 = 90° as are being analyzed here. 




Another view of the evolution of the backwards propagated trajectories from Figs. [25] -[27] are shown 
in Fig. 28 In this figure, we see the absolute value of <f>2 plotted for every trajectory for each crossing. 
The lower plot shows how many trajectories still exist. The main takeaway from this plot is that most the 
trajectories that we propagated backwards don't stay in a librational state for very long; of the 1681 initial 
trajectories, only about 100 cross 4>2 = more than 100 times. Only 2 trajectories last longer than about 
400 crossings. Interestingly, one trajectory lasts for a very long time - approximately 22,000 crossings! This 



tells us that a large portion of the red points seen in Fig. 27 can be attributed to one very rare case (this 



trajectory is shown in Fig. 29 1. The other thing we find from this plot is that, unlike in the u ~ 1 case, 
there isn't a clear value of 4>2 that indicates if circulation will occur or not. This is due to the coupling 
between the libration and orbital states. 

Finally, it is interesting to see how things change when the energy level is increased. Fig. [30] shows the 
02 = ellipsoid for three increasingly higher energy levels. These views show that as the energy is increased 
there are two main impacts. First, the trajectories evolve to much larger areas of ellipsoid, meaning there are 
smaller sets of initial conditions at 02 = that correspond to bounded trajectories. Second, the unbounded 
trajectories reside for much shorter timespans around 02 = 0. In Fig. ]3T] we show the crossing history, 
which for that case has all trajectories circulating again after less than 30 crossings. We also show a stable 
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Fig. 28 Absolute values of 02 for each crossing (those with a negative (f>2 are plotted in blue), along with the number of 
trajectories which are bound at that crossing (out of the initial 1681). The left plot is zoomed in on the first 1000 crossings 
of the right plot; note that after approximately 700 crossings only one sample remains bound. It is clear from the right 
plot that this one particular long-bound case reaches much lower (j>2 values during some crossings than any of the other 
samples. 
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Fig. 29 The trajectory of the long lived case from Fig. |28| in r ~ 4>2 space. The initial trajectory from 02 = 90° — > is 
shown in black, and the final portion of the trajectory is shown in red, before it becomes unbound, in this case toward 
02 = -90°. 



trajectory, proving that although the ZVC is significantly opened at this energy level, there are still bounded 
trajectories. 

7 Conclusion 

This paper explored the dynamic system of a triaxial ellipsoid satellite in orbit in the equatorial plane 
of an oblate body. This system reduces to a 2 degree-of-freedom system in the radial separation and 
the libration angle of the ellipsoid. The reduced system was then analyzed with the goal of determining 
limits on the dynamic configurations for which the librational motion is bounded to less than ±90°. The 
relative equilibrium points were found in the ellipsoid-fixed frame, and the stability of these points was 
discussed. The conservation of energy and angular momentum in the system was exploited to write an 
expression for zero- velocity curves. These curves are used to determine a sufficiency condition on when 
the librational motion is bounded. It is shown that bounded trajectories exist beyond these sufficiency 
conditions, including families of periodic orbits. We addressed the conditions for unbounded motion for 
all systems. In particular, when the ellipsoid becomes of very small mass compared to the oblate body, 
analytical relationships were derived that determine the maximum libration angle for any orbit eccentricity. 
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Fig. 31 An example of a bounded trajectory and the crossing history for the case when E = E + 2.55E. The trajectory 
had initial conditions of r ~ 10, 02 = 0, r ~ 0.02 and ij>2 = 0.5 deg/s and was integrated for 1000 librational periods. 



Future work includes extending this approach to two ehipsoid systems, non-conservative systems, and 
coupled out-of-plane motion. 



Appendix 

The mathematical description of the system studied in this paper was determined as a simplified form of 
the system studied by Scheeres [19] ■ The simplifications made to obtain the results used in this paper are 
derived here in order to make the connection to the previous work explicit. 



Equations of Motion 

The Lagrangian, L = T — V for this system is, 

L = llu^ [e + ,j>,y + ^hj^ + ^ur^ + I (j2^ + ur^) +12^20 - V{r,<j,2) (94) 
Using Lagrange's equations with out any external forces, 

d_ rdL\ _ dL 

dt \dqj dq, ^^^^ 

the equations of motion for this system with the coordinates r, 6^, <j!>i, and <j!>2 were found in Scheeres [19] 
and are rewritten in normalized units to be, 

r = Q r — (96) 

V or 
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However, in the case of an oblate body the potential is no longer a function of 4>i so that 



(100) 



and the equations of motion for the angles become, 



1 dV re 



1 + =r- 



^\ I dV 



I2 J vr"^ 902 



J_dV_ _ ^r6 



+ 2 — 



The partials of the potential for an oblate primary are, 



^ = ^/i + A 

dr r'^ 1 2r^ 



(^1. - Is) - \l2^ - \l2, +l2^+\ - /2J COS 2,^2 



[l2y - l2^) sm(202) 



d(l>2 2 r3 



(101) 
(102) 
(103) 

(104) 
(105) 



Integrals of Motion 

In the current problem with an oblate primary, we have three integrals of motion. The first is the total 
energy of the system, which is shown to be the Jacobi integral of this system since it is time invariant. 



h = ^.^-L = T + V 
oq 



(106) 



The second integral of motion is the total angular momentum of the system. This is found because the 
coordinate 6 is ignorable, meaning that d/dt{dL/d6) = 0, so the integral is written as. 



Kt. 



dL 

de 



Mi- 



(107) 



/ze + /2>2 + ^/l.01 
M2 



The third integral is found by combining Eqs. ( 101 1 and ( 103 ), so we find that 



+ 6» = 



(108) 



Therefore the inertial angular velocity of the primary, 9i, is an integral of motion. This implies that 
the terms in the kinetic energy and angular momentum expressions which depend only on 61 are also 
conserved. This fact makes intuitive sense as a primary that is symmetric about the spin axis can't have 
any gravitational torques exerted on it from the secondary since the center of mass and the center of gravity 
(in the secondary's gravity field) are at the same location in the primary. 
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Dynamics Matrix Partials 



The partials are 



dr 
dr 



di" 
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V drd4>2 
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And the second partials of the potential are given by, 
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(110) 
(111) 
(112) 
(113) 
(114) 
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Relative Equilibria Locations 

Following Scheeres [19 , we find the equilibrium points by searching for places where the variations in 
energy are stationary at a constant value of angular momentum. In other words, we find when the following 
conditions hold: 

(119) 
(120) 
(121) 
(122) 

(123) 
(124) 
(125) 



at a given value of K, as seen in Eq. ( 15 ) 



The following relationships are found. 
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902 2 r3 

which when combined with the stationarity conditions imply, respectively, that f = 0, 02 = 0, and 02 = 0, 
±7r/2, or TV. Using these conditions, we can evaluate the partial with respect to r to be, 



dE_ 
dr 



1 + ^ 



2r2 



(126) 
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Relative Equilibria Stability Partials 



Recall from Sectionjsj the equilibrium points must occur at f = 0, <^2 = 0, and <^2 = 0, ±7!"/2, or n. Therefore 
the partials for the dynamic matrix are greatly reduced to become, 
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Where the equilibrium point distance is indicated by reg. Note that from Eq. (1051, we see that 

^] =0 



(127) 

(128) 
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(130) 
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(134) 



The second partials of the potential evaluated at the equilibrium points are, 



Qj.2 



1+ — 

' eq 

drd<j)2 
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(135) 
(136) 
(137) 



where was defined in Eq. (21 1. The plus terms correspond to the case where <^2 = or tt, and the minus 
terms correspond the the other cases where 4>2 = 7r/2 or 37r/2. 
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